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Abstract 

This paper constructs a doubly robust estimator for continuous dose-response estima¬ 
tion. An outcome regression model is augmented with a set of inverse generalized propen¬ 
sity score covariates to correct for potential misspecification bias. From the augmented 
model we can obtain consistent estimates of mean average potential outcomes for distinct 
strata of the treatment. A polynomial regression is then fitted to these point estimates to 
derive a Taylor approximation to the continuous dose-response function. The bootstrap is 
used for variance estimation. Analytical results and simulations show that our approach 
can provide a good approximation to linear or nonlinear dose-response functions under 
various sources of misspecification of the outcome regression or propensity score models. 
Efficiency in finite samples is good relative to minimum variance consistent estimators. 
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1 Introduction 


The typical set up in causal inference problems is one in which the data available for estimation 
are realisations of a random vector, Zj = {Yi,Di,Xi), i = 1, ...,n, where for the Tth unit of 
observation Yi denotes a response, Di the treatment received, and Xi a vector of pre-treatment 
covariates. The objective is to estimate an average treatment effect (ATE), or in other words, 
the difference in response that would occur under different treatment assignments averaged 
over the population. However, since the treatment is usually not assigned randomly, simple 
comparisons of mean responses across different treatment groups will not in general reveal a 
‘causal’ effect due to confounding. 


Consistent causal estimates of ATEs can, however, be obtained if the vector of covariates Xi is 
sufficient to ensure con ditional independence of pot ential outcomes and treatment assignment. 
Using the notation of Tsiatis and Davidian ( 2007l f. we write joint densities of the observed 
data in the form 

fz{z) = fY\D,x{yi\di,Xi)fr)\x{di\xi)fx{xi)- 


Well known causal ATE estimators can be derived via an outcome regression (OR) model 
for E(l)|Dj, Aj), the mean of the conditional density of the response given the covariates and 









treatment status; or via a propensity score (PS) model, 7r(Z)j|Xj), for the treatment assignment 
(or exposure) mechanism. 


Doubly robust (DR) estimators combine the OR and PS models, usually by weighting or 
augmenting the OR model with covariates derived from the inverse propensity score, to 
obtain estimates of ATEs that are consistent and asymptotically normal when either the 


context of binary treatments (e.g. I 

lobins 2000l. Robins et al. 

2000. Robins and Rotnitzkv 

2001, 

van der I^aan and Robins 200.11 

lainceford and DavidianI 

2004. Banff and Robins! 

2005, 

Kane and Schafei]|2007|l. 


In thi s paper we extend the binary DR augmented regression approach of IScharfstein et al. 


(Il999ll ;o derive a Taylor approximation for continuous dose-response functions. We present 
theoretical results which explain the DR properties of our approach and simulation results 
which show that the estimator can be used to recover a good approximation to linear or 
nonlinear dose-response functions. 


The paper is structured as follows. Section two explains some general principles of PS based 
continuous dose-response estimation. Section three presents our DR model for continuous 
treatments and demonstrates the consistency properties of this estimator. Simulation results 
are presented in section four. Conclusions are drawn in the final section. 


2 Potential outcomes and dose-response estimation for contin¬ 
uous treatments 


Throughout the paper we use the convention of upper-case to denote a random variable and 
lower case for its observed value. In the single time-point setting, let d € D T M, denote a 
given value (dose) of treatment and let d* denote the value of the treatment actually assigned 
to unit i. For each unit i let there be a potential outcome l)(d) associated with each dose, 
such that Ti = {^i(d) : d G V} denotes the full set of potential outcomes. In the data we 
observe random variables describing the actual dose received Di , the outcome associated with 
that dose Yi{Di), and a set of pre-treatment covariates Xj. We define the indicator for receipt 
of dose d as 


Id{Di) 


1 if A = d 

0 if A / d ’ 


and similarly the indicator for receipt of the dose actually received is A (A)- We assume 
that the relationship between actual and potential outcomes obeys the stable unit treatment 
value assumption (SUTVA) such that 1) = /rf(A)Pi(d) for all d G D, for all Tj(d) G A, and 
for i = 1 ,..., n. 


Our starting point for the analysis of continuous treatments within a pot ential outcome s 
framework is the assumption of weak conditional independence introduced by ImbensI ( 2001)1 1. 
which requires that Ti(d) X /^(A)|Xj for all d G D. If weak conditional independence holds 
we can obtain causal estimates from the observed data without knowledge of the full range of 
potential outcomes because 


E[Yiid)\Xi]=E[Yi{d)MDi),Xi]=lK[YMDi),Xi\ 


( 1 ) 
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and taking the expectation over covariates Xi gives the Average Potential Outcome (APO) 
at dose d, fi{d) = E \Yi{d)], commonly referred to as the dose-response function. 

An OR based approach could be used to estimate ([I]) by specifying a mean response model 
E(li|Z)j, Aj) = /3)}, for some regression function m() with known link function 

^ and unknown parameter vector /3. If the assumed model provides a correct specification of 
the true mean response then in principle a consistent estimate of the dose-response at dose d 
can then be obtained using 


n 

J^OR{d) = n~^ ^ |m(d, A*; 


2 = 1 


( 2 ) 


that is, the average across all data of the predicted values evaluated at Di = d 


As in the case of binary treatments, however, estimation of continuous dose-response functions 
is often made more tractable by working with a scalar PS rather than the potentially high di¬ 
mensi onal covariate vector. For continuous treatments, Imbens ( 20001 1 and Hirano and Imbens 
( 2004l i introduce the generalized propensity score (GPS). Let r{d,Xi) = fD\x{d\xi) denote 
the conditional density function for receiving a particular dose of the treatment given pre¬ 
treatment variables Aj = Xi. The observed GPS, which we denote Ri = r{Di,Xi), is a random 
variable comprising values from this conditional density evaluated at the level of treatment 
actually received (i.e. Di) given Aj = x. In addition, we also define the family of random 
variables indexed by d, Rd,i = r{d,Xi), as values from the conditional density of receiving a 
particular dose d given Aj. Clearly when Di = d, Ri = Rd,i- 


Hirano and ImbensI ( 20041 1 show that the GPS provides a bias removal strategy in the context 
of continuous treatments such that 


fi{d)=E[Y,id)]=Ex[E(Y,{d)\X)]=ExmY,{d)\r{d,Xi)]]=ExmY,\Id{D,),r{d,Xi)]]. (3) 


Dose response estimation can therefore proceed via OR or GPS based models. The consistency 
of the estimator (j2]) relies on correct specification of an OR model while that of ([3]) relies on 
correct specification of the GPS. In the next section, we develop a DR estimator for continuous 
and potential ly nonlinear dose- r espon se functions by extending the DR approach for binary 


treatments of Scharfstein et al. ( 1999l l 


3 Doubly robust estimation for continuous treatments 


Dose-response estimation is challenging because the relationship may be nonlinear and the 
bias induced by sources of misspecification non-constant across doses. In constructing an 
appropriate dose-response estimator, flexibility is required to address the inherent trade off 
between fidelity to these nonlinearities and consistency of the estimated causal quantities. 

Our proposed approach involves two main steps. First we specify an augmented OR model 
to obtain consistent DR estimates of mean APOs for strata of the treatment. We then use 
polynomial regression to recover point estimates of the dose-response curve with standard 
errors calculated via the bootstrap. 
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3.1 Estimating equations for DR mean APO estimation 


In the first step, our estimand of interest is the mean APO for strata of the treatment. We 
proceed as follows: 

(i) We define Q, q = (1, strata over the range of d and use T>q to denote treatment 

stratum q {'Dq C P C M). We index distinct treatment levels within each stratum by 
j = (1,J), such that individual elements of Vq are denoted dqj. Our objective is to 
estimate 

fl{Vq)=K,[E{Y,{dq,)}]. (4) 

(ii) Given some assumed density, /£)|x(d|a^i, a), with parameters vector a, and for some 
small 6 defined around the dose of interest, we define a probability representation of the 
GPS as 


7r(d|Aj; a) = Pr(d — 6 < Di < d + 6\Xi,a) 

rd-\-S 

fD\x{t\xi,a)dt 


f 


' d—5 

= E[Id{Di)\Xi = x], 


(5) 


which we refer to as the probabilistic generalized propensity score (PGPS). 

(iii) Specifying an appropriate regression model for /£)|x(d|a:j, a), we use observed doses di 
and observed covariates Xi to obtain consistent estimates of a and calculate observed 
PGPSs: ^{Di\Xi-,a). 

(iv) Denoting membership of treatment stratum q using the indicator Iq{Di), we specify an 
augmented OR (AOR) model as 


e{Di,Xi,Ki{Di,Xi)-^} = ^ ^ 


Q 


m(A,Ai;/3) + ^ 


9=1 


^qIq{Di) 

^{Di\Xi-,a) 




( 6 ) 

say, where ip = (ipi, ...,ipQ) is a Q dimensional parameter vector for the inverse PGPS 
covariates, ^ = (/3, ip), and Ki{Di,Xi) is an (n x Q) matrix each column of which contains 
a covariate for treatment stratum q 


4(A) 

^{Di\Xf,a)' 

We obtain estimates of ^ = {/3,p) as solutions to estimating equations of the form 


E l de {di, Xi,Ki{di, Xi); C} 

i-i ^ 


[Vi 


e{di,Xi,Ki{di,Xi;^)}] = 0 


(7) 


where (/> is a working variance matrix for Var[li| Aj A,, Kj(A) Aj)]. 

(v) For each stratum, we then calculate the means of the predicted values of the AOR model 
evaluated at each level of 'Dq. Our estimator is the average of these means 


d-DR (Dq) 



Vq \ 

'^idqjlXi, Cx) j 


( 8 ) 
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(vi) Finally, we slide the partition to define new strata and repeat steps (i) to (v) thus 
increasing coverage of mean APOs on the dose-response function. 

For valid inference, the calculations must be based on a defined region of dose, CCD, 
within which there is common support by treatment status in the covariate distributions. The 
common support requirement is met when for any subset of C, say ^ C C, Pr((i G A\Xi = 
x) > 0 for all d and x and A C C. 

The estimator is doubly robust in the sense that it yields a consistent estimate of the mean 
APO in each stratum if one of two parametric restrictions hold: 1) the OR model Di'A)} 

provides an consistent estimate of the true conditional expectation E(yi|T)j, Aj), or; 2) 7f(d| Aj; a) 
provides a consistent estimator of the true POPS E[/rf(T)j)|Aj = x\. These restrictions, and 
their implications for mean APO estimation, are described in the following two theorems on 
redundant conditioning and bias correction. 

Theorem 1. (Effects of redundant conditioning on APO estimation). If the OR model 
'h“^{m(Aj, D,;/3)} provides a consistent estimate o/E[li|Dj, Aj], then the mean of the pre¬ 
dicted values from this model evaluated at Di = d will provide a consistent estimate of¥.\Yi{d)] 
provided conditional independence between treatment status and response holds. The mean pre¬ 
dicted values from the AOR model, evaluated at Di = d, will also produce a consistent estimate 
ofK[Yi{d)] under conditional independence, but the AOR model will be less efficient. 

Proof. The true conditional expectation to be estimated is 'K\Yi\Di, Xf, but instead we esti¬ 
mate a model for K Aj,/?j(Dj, Aj)]. If we have conditional independence given covariates 

Xi, and if the SUTVA holds, then 

E[Yi{d)\Xi,Ki{Di,Xi)]=K[Yi\d,Xi,KiiDi,Xi)], 

regardless of redundant conditioning on'ki{Di,Xi). The dose-response can then be obtained 
by taking expectations 

E[YM)\=^x,,n,{D,,xM^M,Xi,Ki{Di,Xi)]]=ExAny^\d,Xi]\. 

However, the AOR model will have higher variance than the OR model because 

Var (E [y, id)\Xi , (A, X*)]) =Var (E [E [A, (d) | A, % {Di , A)] (ki (A, A)]) 

+ E [Var (E [V(d)|A, k,(A, A)] |k,(A, A))] 

=Var (E [yi(d)|A]) +E[Var (E A(d)|A, Ki(A, A)] \ki{D„Xf))] 

and since E [Var (E [yj((i)IAj, Kj(A) A)] |Ki(A) A))] > 0 , tden Var (E[yj(d)|Aj, Ki(A, A))] > 

Var(E[yi(d)|A]). 

Theorem 2. (The bias correction property). If the estimated propensity score Tr{d\Xi;a) 
provides a consistent estimate of the true propensity score 7r(d|Aj;a) = E[/rf(A)|A = x\, 
then the mean predicted values from the AOR model, averaged within some stratum of the 
treatment, will provide consistent estimates of the mean APO for that stratum, fj,{'Dq), under 
mis specification of the OR model. 


Proof. Parameter estimates (p are obtained from the AOR model as solution to estimating 
equations of the form 


n 




4(A) 


^ 7fi(A|A;S) 


Q 

Vi - { m {di,Xi; fi)+'^(pq 

g=l 


Iq{di 


TTi {di \Xi , o) 


= 0 , 
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which converge in probability to 


E 


1 




^i{Di\Xi-,a) 

Writing the AOR model as '1'“^ {•} and taking expectations over Xi and Di gives 


Ex, 


Di 


E /d.(A 


1 


• {•}] 


^i{Di\Xi-,a 

From weak conditional independence and the SUTVA 

1 


Iq{Di) — l,Xi,Di 


^Xi,Di 


E[Irf,(A)|X,]-4(A)-E 


[y,(A) 


^ TTj (^Di \Xi , (X 

and since Xj is a function of Xi we use the pull through property to give 

1 


Iq{Di) — 1, Xi 


^Xi,Di 

= EXi,D. 


e[A(A)|A]- 

7ri(A|A;a) 


^i(A|A;a) 

1 

7ri{Di\Xi;a) 


•4(A)-Efyi(A)-^"'{-} 


Ig{Di) = l,X., 

IqiDi) = 1, Xi 


= ED; 


E y,(A)-^~U-} 


IqiDi) = 1 


= 0 . 


We can write this equation for each stratum q in the form 


PL {Vq) = Ej [Ej {Yi{dqj)}] = Ej 


E ^ < m {dqj, Xi]/3) + 


Fq 


x{dqj \Xi^ OL 


Thus, averaging the mean values of the AOR model over dqj G Vq provides an unbiased 
estimate of p{Vq) regardless of whether '1'“^ |m ^A, A;/3^| is unbiased for E[Yi\Di, Xi], 


The estimating equations of the DR model are implied by a number of standard estimators. 
For instance, Maximum Likelihood Estimation (MLE), Maximum Quasi-Likelihood (MQL), 
Restricted MLE (REML) for linear mixed models (LMMs), and Penalized Quasi-Likelihood 
(PQL) for generalized linear mixed models (GLMMs) 


3.2 Doubly robust dose-response estimation via Taylor series expansion 

The approach outlined above yields consistent estimates of the mean APO for each stratum. 
To obtain point estimates on the continuous dose-response function we apply a Taylor series 
expansion via a polynomial regression model. 

The true dose response function, fj,{d) = E[yj((i)] = f{d]9) is unknown, but we assume that 
it can be well approximated by an Mth degree polynomial function 

M 

E[y*(d)] ~ 

m=0 
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The mean APO for stratum q relates to the polynomial dose-response in the following way 

M 

m=0 


with sample analogue 

M J im 

(9) 

m=0 j=l 

where J is the number of units in stratum q. To evaluate this expression we need estimates of 
6m- We obtain these by regressing the mean APO estimates for each stratum, JliVq), on M — 1 
covariates Ylj = (1) Use of OLS will yield best linear unbiased estimates of 

9m under the Gauss Markov conditions. We use then use the OLS estimates of dm to recover 
a polynomial approximation to points on the dose response function via 


M 

m=0 


Variances of the resulting dose-response point estimates can be obtained by applying a single 
bootstrap resampling scheme over the estimation of each model component (i.e. PGPS, AOR, 
polynomial regression). 


3.3 Comparison with existing doubly robust approaches 


To our knowle dge there are currently only tw o papers that have constructed DR estimators 
using the GPS (jvan der Laan and Robinsll2nn,'ll . construct a DR estimator for continuous treat¬ 
ments but in the context of marginal structural models and not for point exposure problems). 


Flores and Mitnik ( 20091 1 combine OR with inverse GPS weighting to obtain a DR weighted 
regression estimator for multi-valued, but not continuous treatments, in which doses take 
discrete values d = 1, ...k. They do this by weighting the OR model 


k 

E[y,| A, = ii(A = j) + 

j 


with weights Wi = ^JXjRi. To calculate a weighted regression (WR) estimator at dose d they 
use 

Ya= 1 + ■ Wi 

Kd)wRi = --• ( 10 ) 

Ei=i Wi 

They do not provide a formal proof for the DR property but make two arguments for il¬ 
lustration. First, that weighting does not affect the consistency properties of the outcome 
regression; and second, that weighting by a correctly specified GPS gives a balanced sample 
analogous to adjusting for observed characteristics with experimental data. 


construct a similar DR estimator for continuous treatments by combining 
model through inverse probability weighting. Their weighted regression 

estimator is 

1 ^ — 

k‘{d)wR2 = ^^rn{d,Xi,P*), ( 11 ) 

i=l 


Zhang et al.l ((2013) 


and OR and a PS 
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where (3* is the estimate for the weighted regression model m{Di,Xi,l3*) with weights W(Di)/Ri 
where W{Di) > 0 is chosen to stabilize the weights when the GPS takes very small values. In 
their application they chose W{Di) to be the marginal density of Di under a normal model. 

The key difference between our augmented GPS approach, and the GPS weighting approaches, 
is that by including multiple inverse GPS covariates for different strata we are able to induce 
local rather than global bias correction. This follows from theorem 2 above where we show that 
the bias correction is specific to the mean APO for each stratum. If instead we were to employ 
a single weighting scheme, or include a single inverse GPS covariate, the bias correction would 
be tailored to a mean APO over the whole sample. For non-linear dose response-functions, 
local rather than global approximations allow to us to correct for misspecification due to an 
inaccurate assumed functional form in addition to effects of confounding. We illustrate these 
properties in the simulations presented in the next section. 


4 Simulations 


The model structure analysed in our simulations corresponds to the following data generating 
process. 


A continuous treatment is assigned as a function of covariates Xi, X 2 , and Z. The Gaussian 
dose-response function is quadratic in the treatment with confounding from Xi and A 2 . 

Ai, A 2 ~ = 4, /iX2 = 8, cr^i = 1, o-\^ =2,p = -0.5) 

Z ~ AA(10,4) 

D ~ A/'(2.0 -h 0.5Ai -h O. 25 A 2 -h Z, erf,) 

Y ~ A/'(1.0 -h 4.0D - 0.125T)2 -h 0.5Ai 2 A 2 - 0.5A|, cj^) 

The correct OR model is 


E[Y\D, Ai, A 2 ] = /3o + /3iT> + /32F»' + + 1 ^ 4 X 2 + ^^Xl 

and the correct PGPS model is 


1 1 f 

= / - exp 

JD-5 D2TTa‘jj ^ 


1 


2al 


{t - pof ] dt 


where parameters pD and a?) are estimated from the correct regression model 


E[Z)| Ai, A 2 ] — do + aiAi -|- a2A2. 


4.1 Simulations for strata mean APO estimation 


The following models are tested: 


1. piT>k)oRi - based on the correctly specihed OR model with predicted values averaged 
over j = (1,..., J) treatment levels within each stratum 


iioRi{T3>k) = d ^ 


n 


j=l L i=l 


n 







2. 'fl(T>k)oR 2 - same estimator as 'fl{T>k)oRi but based on an incorrectly specified OR model, 
with Xi assumed as sole confounder. 

3. 'fi{Vk)DRi - a DR estimator as described in equation ([8|) above based on an incorrectly 
specified OR model {Xi as sole confounder) augmented with correctly estimated inverse 
POPS (vf^ ) covariates for defined strata of the treatment. 

4. 'j2{'Df;)DR2 - a DR estimator based on a correctly specified OR model augmented with in¬ 
correctly estimated inverse POPS covariates with Xi assumed as sole confounder, 

for dehned strata of the treatment. 

5. 'fl(T>k)DR 3 - a DR estimator based on an incorrectly specihed OR model augmented with 
incorrectly estimated inverse POPS covariates. 

6. 'fi{Vk)DR 4 : - a DR estimator based on an correctly specified OR model augmented with 
correctly estimated inverse POPS covariates. 


The simulations are based on 1000 runs on generated datasets of size 10,000. The mean of d is 
16 and the range approximately 1 to 30. Estimates are presented for the following treatment 
strata: (10.5,12.5], (12.5,14.5], (14.5,16.5], (16.5,18.5], (18.5,20.5]. To calculate the PGPSs we 
set 5 = 0.5. Table [1] shows our simulation results. The following results are reported: average 
estimates (Av Est), average estimated variances (Av Est Var), empirical variances (Emp Var) 
calculated from the bootstrap replications, mean squared error (MSE), and coverage based on 
normal bootstrap 95% confidence intervals. 

The correctly specified OR model, ORl, provides unbiased estimates of points on the quadratic 
dose-response and consequently shows good results for all treatment intervals with relativity 
low variance and good coverage. The incorrectly specihed OR model, OR2, fails to hnd 
the quadratic relationship between outcome and treatment levels, instead indicating a linear 
decreasing effect with poor coverage and large bias and MSE. Augmenting the incorrect OR 
model with the inverse POPS covariates substantially corrects for bias from confounding 
and from functional misspecihcation of the treatment covariate. The DRl estimator hnds 
a quadratic relationship between outcome and treatment and performs well in terms of bias, 
MSE and coverage, although as we would expect, the variances are larger than for the correctly 
specihed OR model (which has minimum variance amongst all consistent estimators). The 
DR2 model shows that the inclusion of addition irrelevant covariates in the AOR model does 
not induce bias, but as illustrated in theorem 1, it does increase the variance relative to the 
correct OR model. If the POPS model and the OR model are wrong there is no DR property 
to be had and this is demonstrated in the DR3 results which indicate large bias and variance 
in estimation of the dose-response and poor coverage. Einally, if both the PS and the OR 
model are correctly specihed the doubly robust estimator (DR4) provides unbiased estimates 
but with variance larger than would be achieved via a correct OR model alone. Thus, in 
summary, the hnite sample properties of our DR approach indicate consistent estimation with 
relatively good performance in terms of efficiency. 


4.2 Simulation of existing doubly robust weighting approaches 

Table [2] shows comparative results obtained by htting the models of equation [10] (WRl) and 
equation [11] (WR2). The hrst thing worth noting is that when the OR model is correctly 
specihed (DR2 models) both WR models provide consistent estimates, and this is because 
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Table 1: Simulation results for Gaussian dose-response GLM with quadratic treat- 
ment effect._ 



(10.5,12.5] 

treatment intervals 
(12.5,14.5] (14.5,16.5] 

(16.5,18.5] 

(18.5,20.5] 

Truth 


15.426 

17.176 

17.926 

17.676 

16.426 

ORl 

Av Est 

15.425 

17.176 

17.926 

17.677 

16.428 


Av Est Var 

0.012 

0.010 

0.010 

0.010 

0.011 


Emp Var 

0.011 

0.010 

0.010 

0.010 

0.011 


MSE 

0.012 

0.010 

0.010 

0.010 

0.011 


Coverage 

95.400 

95.200 

94.800 

94.500 

95.200 

OR2 

Av Est 

16.944 

16.628 

16.311 

15.994 

15.677 


Av Est Var 

0.025 

0.015 

0.010 

0.012 

0.019 


Emp Var 

0.016 

0.007 

0.004 

0.008 

0.019 


MSE 

2.320 

0.308 

2.614 

2.838 

0.579 


Coverage 

0.000 

0.417 

0.000 

0.000 

0.000 

DRl 

Av Est 

15.430 

17.177 

17.925 

17.681 

16.428 


Av Est Var 

0.072 

0.046 

0.037 

0.038 

0.050 


Emp Var 

0.072 

0.046 

0.037 

0.038 

0.050 


MSE 

0.074 

0.047 

0.040 

0.039 

0.051 


Coverage 

94.500 

94.400 

95.500 

94.400 

94.900 

DR2 

Av Est 

15.428 

17.178 

17.926 

17.679 

16.424 


Av Est Var 

0.024 

0.017 

0.016 

0.016 

0.019 


Emp Var 

0.024 

0.018 

0.014 

0.016 

0.019 


MSE 

0.024 

0.017 

0.014 

0.016 

0.019 


Coverage 

95.833 

94.792 

94.792 

96.250 

94.375 

DR3 

Av Est 

16.162 

17.576 

18.012 

17.439 

15.876 


Av Est Var 

0.072 

0.045 

0.037 

0.040 

0.059 


Emp Var 

0.080 

0.043 

0.037 

0.042 

0.062 


MSE 

0.621 

0.203 

0.044 

0.098 

0.365 


Coverage 

25.833 

51.875 

91.042 

77.708 

38.125 

DR4 

Av Est 

15.430 

17.172 

17.930 

17.676 

16.422 


Av Est Var 

0.024 

0.018 

0.016 

0.016 

0.020 


Emp Var 

0.025 

0.017 

0.016 

0.016 

0.020 


MSE 

0.025 

0.017 

0.016 

0.016 

0.020 


Coverage 

94.100 

93.000 

95.000 

94.800 

94.800 
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theorem 1 above holds. Note, however, that the variance of the WR1-DR2 estimates is rela¬ 
tively large. When the OR model is incorrectly specified however, the WR estimators based 
on a global weighting scheme produce biased estimates of the mean APO for each stratum. 
Due to the inclusion of dummies in WRl-DRl model, a quadratic shape is implied in the 
strata estimates, but the global bias correction induced by weighting does not adequately 
address local misspecification. With no dummies, as in WR2-DR1, the model fails to find 
the quadratic shape of the dose response and again mean APO estimates are biased. Note 
that excluding dummies from the WRl approach we get estimates that are approximately 
equivalent to those of WR2. 


Table 2: Comparative results based on GPS weighting estimators. 



(10.5,12.5] 

treatment intervals 
(12.5,14.5] (14.5,16.5] 

(16.5,18.5] 

(18.5,20.5] 

Truth 


15.426 

17.176 

17.926 

17.676 

16.426 

WRI-DRI 

Av Est 

15.854 

17.409 

17.965 

17.561 

16.160 


Av Est Var 

0.068 

0.044 

0.037 

0.039 

0.053 


Emp Var 

0.065 

0.044 

0.036 

0.040 

0.056 


MSE 

0.248 

0.098 

0.037 

0.054 

0.127 


Coverage 

63.750 

82.083 

95.625 

91.458 

79.583 

WRI-DR2 

Av Est 

15.420 

17.171 

17.921 

17.665 

16.420 


Av Est Var 

0.023 

0.017 

0.016 

0.016 

0.020 


Emp Var 

0.022 

0.017 

0.015 

0.017 

0.019 


MSE 

0.022 

0.017 

0.015 

0.017 

0.019 


Coverage 

94.792 

95.208 

95.625 

95.625 

95.417 

WR2-DR1 

Av Est 

16.700 

16.525 

16.351 

16.176 

16.002 


Av Est Var 

0.012 

0.009 

0.009 

0.011 

0.015 


Emp Var 

0.012 

0.010 

0.009 

0.011 

0.015 


MSE 

1.635 

0.433 

2.491 

2.261 

0.196 


Coverage 

0.000 

0.000 

0.000 

0.000 

6.875 

WR2-DR2 

Av Est 

15.420 

17.184 

17.941 

17.692 

16.437 


Av Est Var 

0.012 

0.010 

0.010 

0.010 

0.011 


Emp Var 

0.012 

0.011 

0.011 

0.011 

0.012 


MSE 

0.012 

0.011 

0.011 

0.012 

0.012 


Coverage 

94.700 

95.100 

95.000 

95.300 

95.600 


4.3 Testing for joint significance of the AOR inverse GPS covariates 

In contrast to doubly robust weighting approaches, use of an AOR model allows tests to be 
performed on the inclusion of inverse GPS covariates for bias correction. In table E] below 
we present results for Wald Chi-squared tests for the joint significance of the inverse GPS 
covariates in our DR simulations. 

When the OR model is correctly specified (i.e. DR2 &: DR4) the rejection rate of the null 
hypothesis {ip = 0) is approximately 5% given either a correctly or incorrectly specified GPS 
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Table 3: Simulation of Wald Chi-squared tests for joint significance of inverse 
GPS covariates, rejections rates for the null: Hq : ip = 0 (nominal level is 95%) 


Hq rejection rate 

DRl 

100.0% 

DR2 

4.7% 

DR3 

100.0% 

DR4 

4.8% 


model. The size of the test is thus as expected. When the OR model is incorrectly specified 
(i.e. DRl & DR3) the rejection rate of the null is 100%, regardless of whether the GPS model 
is correctly specified or not. Thus, while the test shown here cannot provide inference on how 
well the GPS covariates have helped to reduce bias from confounding, it can help indicate 
when their inclusion offers little beyond an OR specification. This is important, because as 
shown above, inclusion of irrelevant covariates in the AOR model can reduce efficiency. 


4.4 Simulations for dose-response estimation 

We now illustrate how Taylor expansion via polynomial regression can be used to recover an 
approximation to the dose-response curve using our mean APO estimates. For the same set up 
described above we obtain 1000 simulated set of results for different treatment strata, sliding 
the partition of the treatment to increase coverage of mean APOs. We then take the mean 
APO estimates for each strata as a response variable and fit a second order polynomial model 
to the mean dose as described in section 3 above. The parameters of the polynomial model 
are then used to recover an approximation to the dose response. 

The results are illustrated graphically in figure [1] which shows the resulting dose-response 
curve obtained for DRl and OR2 estimates relative to the truth. The figure also shows the 
mean dose-response curve obtained when the OR2 model is fitted using a Generalised Additive 
Model (GAM) (OR2sp), in which case model misspecification arises from confounding rather 
than functional misspecification. 

The nonparametric fit to the DRl estimates recovers a good approximation to the true dose- 
response curve. The OR2 model finds an incorrect shape for the dose-response because it is 
based on biased estimates of the mean APO for each stratum as discussed above. Fitting 
the OR2 model semiparametrically, with a smooth term for treatment d, gives a quadratic 
dose-response curve but which is still biased due to the effect of confounding from omitted 
covariates. This comparison thus demonstrates the ability of the DR model to correct for 
misspecification bias arising from both confounding and an incorrect functional form. 


5 Conclusions 


In this paper we have constructed a doubly robust approach for estimation of dose-response 
functions for continuous treatments. The model is doubly robust in the sense that consistent 
estimates of average potential outcomes can be obtained for defined strata of the treatment 
under misspecification of either the outcome regression or propensity score models. This is 
achieved by inducing local bias correction for misspecification of the outcome regression model 
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Figure 1: Estimated dose-response functions obtained from fitting second order 
polynomial models to the DRl and OR2 mean APO estimates. 


by augmenting it with inverse propensity score covariates corresponding to defined strata of 
the treatment. A polynomial regression is then fitted to these point estimates to derive a 
Taylor approximation to the continuous dose-response function. Standard errors are derived 
by bootstrapping over all model components. 

We have shown that our approach can provide a good approximation to linear or nonlinear 
dose-response curves and is robust to problems of confounding or functional form misspeci- 
fication. Furthermore, simulations show that the efficiency performance in finite samples is 
good relative to minimum variance consistent estimators. There are several issues in the pa¬ 
per that require further attention in future research. In particular, the choice of treatment 
strata for point estimation could be given a more formal basis. Similarly, the choice of delta 
for estimation of propensity score probabilities should be examined in greater depth than has 
been possible here. 
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